/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#define max(x,y) (((x)>(y))? (x):(y))
#define min(x,y) (((x)<(y))? (x):(y))

#include "fdlayout.h"
#include "graph.h"
#include "controller.h"

ForceDirectedLayout::ForceDirectedLayout(const char *name): Layout("ForceDirected", name) {
    setTimestep(0.2);
    setRunning(true);
    timestep = 0.1;
    K = 0.2;  // Elasticity Constant
    T = 0.5; // temperature
    repulsion = 0.02; // Repulsion
    springRest = 0.0;

    setPropertyInt("Iterations", 100);

    setPropertyFloat("Elasticity", 0.2);
    setPropertyFloat("Spring Rest", springRest);
    setPropertyFloat("Repulsion", 0.02);
    setPropertyFloat("Time step", 0.1);
    setPropertyFloat("Temperature", T);
    setPropertyBool("Enabled", false);
    setPropertyFloat("MaxDisplacement", 1.0);
    setPropertyFloat("MaxDistance", 1.0);
    setPropertyStringList("Edge Filter", "Default", "edge");

    timer = createTimer(20, this);
}

void ForceDirectedLayout::setTimestep(float ts) {
    timestep = ts;
}

float ForceDirectedLayout::getEnergy() {
    return totalEnergy;
}

bool ForceDirectedLayout::isRunning() {
    return running;
}

void ForceDirectedLayout::setRunning(bool value) {
    this->running = value;
}

int random(int a, int b) {
    return (int) ((float) b * (float) rand()/(float) RAND_MAX  + a);
}

void ForceDirectedLayout::updateVisible() {
    visible.clear();
    for(int i=0;i<numNodes();i++) {
        if(isNodeinLayout(i)) {
            visible.push_back(i);
        }
        else setForce(i, 0, 0);
    }
}

#define sgn(x) (((x)<0)? (-1.0):(1.0))

void ForceDirectedLayout::computeIteration() {
	// TO DO:
	// Fix the node with largest degree
    assert(graph);
    float maxVelocity = 5.0;
    float maxDisplacement = getPropertyFloat("MaxDisplacement");
    if(!running || timestep<=0) return;

    string edgeFilterProperty = PROPERTY(getPropertyString("Edge Filter"));

    float totalEnergy = 0;
    int N = numNodes();
    for(unsigned int ki=0;ki<visible.size();ki++) {
        int i = visible[ki];
        setForce(i, 0, 0);
    }
	int *degree = new int[graph->nodes.size()];
	for(unsigned int i=0;i<graph->numNodes();i++) degree[i]=0;

    for(unsigned int ki=0;ki<visible.size();ki++) {
        int i = visible[ki];
        for(int k=0;k<numNeighbors(i);k++) {
            int j = getNeighbor(i,k);
			bool useEdge = graph->getAsBoolean(graph->nodes[i]->edges[k], edgeFilterProperty);
			if(useEdge) {
				degree[i]++;
				degree[j]++;
			}
		}
	}
	for(unsigned int ki=0;ki<visible.size();ki++) {
        int i = visible[ki];
		for(int k=0;k<numNeighbors(i);k++) {
            int j = getNeighbor(i,k);
			bool useEdge = graph->getAsBoolean(graph->nodes[i]->edges[k], edgeFilterProperty);
			
            if(isEdgeinLayout(i,j) && useEdge) {
				float fx, fy;
				if(degree[j]>1) {
					getAttractionForces(i,j,fx,fy);
					addForce(i, fx, fy);
				}
				if(degree[i]>1) {
					getAttractionForces(j,i,fx,fy);
					addForce(j, fx, fy);
				}
            }
        }
		for(unsigned int k=0;k<visible.size();k++) {
            int j = visible[k];
            if(i!=j) {
				float fx, fy;
				getRepulsionForces(i,j, fx, fy);
				addForce(i, fx, fy);
            }
        }
    }
    for(int i=0;i<N;i++) {
        if(!isNodeinLayout(i)) continue;
        if(!isNodeUpdatable(i)) continue;
        float fx, fy;
        float vx, vy;

        float newvx, newvy;
        float x, y;
        float newx, newy;

        getForce(i, fx, fy);
        getVelocity(i, vx, vy);
        getLocation(i, x, y);

        newvx = (vx + timestep*fx)*T;
        newvy = (vy + timestep*fy)*T;

        if(fabs(newvx)>maxVelocity) newvx = sgn(newvx)*maxVelocity;
        if(fabs(newvy)>maxVelocity) newvy = sgn(newvy)*maxVelocity;

        float dispx = timestep*newvx;
        float dispy = timestep*newvy;
        if(fabs(dispx)>maxDisplacement) dispx = sgn(dispx)*maxDisplacement;
        if(fabs(dispy)>maxDisplacement) dispy = sgn(dispy)*maxDisplacement;

        newx = x + dispx;
        newy = y + dispy;

        float d = newx*newx + newy*newy;
        float R = 100.5;
        if(d>R*R) {
            newx = R*newx/sqrt(d);
            newy = R*newy/sqrt(d);
        }
        setVelocity(i, newvx, newvy);
        setLocation(i, newx, newy);

        float speed2 = newvx*newvx + newvy*newvy;
        float energy = 1.0*speed2;
        totalEnergy +=energy;
    }


    T = T*0.999; // cooling off
    if(running) {
        //if(timestep<0.1) timestep = 0.1;
    } else {
        if(totalEnergy<1e-4) timestep = 0;
    }
    updated();
	delete degree;
}

int ForceDirectedLayout::numNodes() {
    return graph->nodes.size();
}


void ForceDirectedLayout::getAttractionForces(int i, int j, float &fx, float &fy) { //Hooke Variant
    float l = springRest;  // rest length
    float x1,y1,x2,y2;
    getLocation(i,x1,y1);
    getLocation(j,x2,y2);
    float dx = x1-x2;
    float dy = y1-y2;
	

    float d = (dx*dx+dy*dy);

    float f = -K*(d-l);
    if(d>0) {
        fx = dx/d*f;
        fy = dy/d*f;
    } else {
        fx = 0;
        fy = 0;
    }
}

void ForceDirectedLayout::getRepulsionForces(int i, int j, float &fx, float &fy) {
    float q1 = 1.0; //(float) graph->nodes[i]->getPropertyInt("Degree")/100000.0;
    float q2 = 1.0; //(float) graph->nodes[j]->getPropertyInt("Degree")/100000.0;

    float x1,y1,x2,y2;
    getLocation(i,x1,y1);
    getLocation(j,x2,y2);
    float dx = x1-x2;
    float dy = y1-y2;

    float d2 = (dx*dx+dy*dy);

    float d = sqrt(d2);
    if(d>0) {
		float f = d2>0 && d2<0.5? 0.001*repulsion*q1*q2/(d2+1e-3): 0;
        //float f = repulsion*q1*q2/d2;
        fx = dx/d*f;
        fy = dy/d*f;
    } else {
        fx = fy = 1e-4 + springRest;
    }
}

inline bool ForceDirectedLayout::isNodeinLayout(int i) {
    return graph->nodes[i]->isVisible();
}
bool ForceDirectedLayout::isEdgeinLayout(int i, int j) {
    return isNodeinLayout(i) && isNodeinLayout(j);
}
bool ForceDirectedLayout::isNodeUpdatable(int i) {
    return !graph->nodes[i]->isSelected();
}
int ForceDirectedLayout::numNeighbors(int i) {
    return graph->nodes[i]->edges.size();
}
int ForceDirectedLayout::getNeighbor(int i, int k) {
    return graph->nodes[i]->edges[k]->getNode1()->getIndex();
}

float random(float a, float b) {
    return a + (b-a)*(float) rand()/(float) RAND_MAX;
}
void ForceDirectedLayout::reset() {
    for(unsigned int i=0;i<graph->nodes.size();i++) {
        float x,y;
        x = random(-1,1);
        y = random(-1,1);
        setLocation(i,x,y);
        setVelocity(i,0,0);
    }
}

void ForceDirectedLayout::timerEvent() {
    updateVisible();
    computeIteration();
    controller->updateRender();
}

void ForceDirectedLayout::compute() {
    int numIterations = getPropertyInt("Iterations");
    for(int i=0;i<numIterations;i++) {
        printf("."); fflush(stdout);
        updateVisible();
        computeIteration();
    }
    printf("\n");
}

void ForceDirectedLayout::update(string property) {
    Layout::update(property);
    springRest = getPropertyFloat("Spring Rest");
    repulsion = 0.01*getPropertyFloat("Repulsion");
    K = getPropertyFloat("Elasticity");
    timestep = getPropertyFloat("Time step");
    float temp = getPropertyFloat("Temperature");
    if(temp>T) T = temp;
    bool enabled = getPropertyBool("Enabled");
    assert(timer);
    if(enabled) {
        timer->start();
    } else {
        timer->stop();
    }
}

Layout *createForceDirectedLayout(string name) {
    ForceDirectedLayout *layout = new ForceDirectedLayout(name.c_str());
    return layout;
}
